!
!  FLAT     Flat field correction utility.
!  Copyright (C) 1997-2005  Filip Hroch, Masaryk University, Brno, CZ
!
!  $Id$

Subroutine Flat( sci, flt, out, nflat)

  implicit none

  character(len=*) :: sci, flt, out
  integer :: nflat

! Parameters

! No. of image dimensions
  Integer, parameter :: DIM = 2

! The flatbat assumes that output image has in intensity
! within range minvalue < x < maxvalue.
! The minvalue always is 0 for measured intensities (?)
! The maxvalue is determined as (2**bitpix - 1) or by SATURATE
! keyword in image headres (in specified order).
! The specific characteristic of dark images may require some modification
! for real-type intensity of image.
  Real, parameter :: minvalue = 0.0
  Real            :: maxvalue = 65535.0

! FITS variables for flatimage
  Integer :: fstat,fBlockSize,fNaxis,fNaxes(DIM),fPCount,fGCount,fBitPix
  Logical :: fSimple,fExtend
  Real :: fscale, fzero, bscale, bzero
  Real, allocatable, save :: FLATDATA(:,:),fbuf(:)

! FITS variables for scientific exposition
  Integer :: stat,sBlockSize,sNaxis,sNaxes(DIM),sPCount,sGCount,BitPix
  Logical :: sSimple,sExtend
  Real :: scale, zero
  Real, allocatable, save :: SCIDATA(:,:)

! FITS variables for output image
  Integer :: ostat
  Real, allocatable, save :: OUTDATA(:,:)

  Logical :: exist, anyf
  Integer :: i, j, nmean
  Real, save :: xmean, xsig
  Real :: x
  Character(len=256) ::  buf, scalecom, zerocom

!-----

  fstat = 0
  stat = 0
  ostat = 0
  
      if ( .not. allocated(FLATDATA) ) then  
 
! flatframe initialisation
        Call ftopen(20,flt,0,fBlockSize,fstat)
        Call ftghpr(20,DIM,fSimple,fBitPix,fNaxis,fNaxes,fPCount,fGCount,fExtend,fstat)
        if( fstat /= 0 ) goto 666
        Call ftgkye(20,'BSCALE',fscale,buf,fstat)
        Call ftgkye(20,'BZERO' ,fzero ,buf,fstat)
!        BSCALE, BZERO not presented
        if( fstat /= 0 ) then
           fstat = 0
           bscale = 1.0
           bzero = 0.0
        endif
!        parameters
        Call ftgkye(20,'SKYMOD',xmean,buf,fstat)
        Call ftgkye(20,'SKYSIG',xsig ,buf,fstat)
        if( fstat /= 0 ) then
           fstat = 0
           xmean = 0.0
           xsig = -1.0
        endif

        Allocate(FLATDATA(fNaxes(1),fNaxes(2)))
        Call ftg2de(20,1,minvalue,fNaxes(1),fNaxes(1),fNaxes(2),FLATDATA,anyf,fstat)
        Call ftclos(20,fstat)
        if( fstat /= 0 ) goto 666
        if( anyf ) write(*,*) "Warning: undefined pixels."

! median computation
        allocate(fbuf(fNaxes(1)*fNaxes(2)))
        nmean = 0
        do i = 1, fNaxes(1)
          do j = 1, fNaxes(2)
             nmean = nmean + 1
             fbuf(nmean) = flatdata(i,j)
          enddo
        enddo
        Call RobustMean(fbuf,nmean,xmean,xsig)
        deallocate(fbuf)
        if( xmean /= 0.0 )then
          FLATDATA = FLATDATA / xmean
        else
          FLATDATA = huge(0.0)
        endif

      endif

!      write(*,*) & 
!        '(',trim(sci),') / (',trim(flt),'/',xmean,') = (',trim(out),')' 

! scientific exposure initialisation
      Call ftopen(25,sci,0,sBlockSize,stat)
      Call ftghpr(25,DIM,sSimple,BitPix,sNaxis,sNaxes,sPCount,sGCount,sExtend,stat)
      if( stat /= 0 ) goto 666
      Call ftgkye(25,'BSCALE',scale,scalecom,stat)
      Call ftgkye(25,'BZERO' ,zero ,zerocom,stat)
!        BSCALE, BZERO not presented
        if( stat /= 0 ) then
           stat = 0
           bscale = 1.0
           bzero = 0.0
        endif
      if( .not. allocated(SCIDATA) ) Allocate( SCIDATA(sNaxes(1),sNaxes(2)))
      Call ftg2de(25,1,minvalue,sNaxes(1),sNaxes(1),sNaxes(2),SCIDATA,anyf,stat)
      if( stat /= 0 ) goto 666

! a range of intensities of the image
     maxvalue = 2.0**BitPix - 1.0
     call ftgkye(25,'SATURATE',x,buf,stat)
     if( stat == 0 ) maxvalue = x
     stat = 0

! output initialisation - parameters by scientific image
      Inquire(file=out,exist=exist,err=666)
      if( exist )then
        open(26,file=out,err=666)
        close(26,status='delete',err=666)
      endif

      if( .not. allocated(OUTDATA) ) Allocate( OUTDATA(sNaxes(1),sNaxes(2)))
      Call ftinit(26,out,sBlockSize,ostat)
! header keywords copying - by scientific image
      Call ftcopy(25,26,1,ostat)
      Call ftphis(26,'Flatframe corrected.',ostat)
      if( ostat /= 0 ) goto 666

! flat divide
      where( FLATDATA /= 0.0 )
          OUTDATA = SCIDATA / FLATDATA
      elsewhere
          OUTDATA = maxvalue
      end where

! range checking
!      where( OUTDATA < minvalue )
!        OUTDATA = minvalue
!      endwhere
!      where( OUTDATA > maxvalue )
!        OUTDATA = maxvalue
!      endwhere

! corrected image output
      call ftp2de(26,1,sNaxes(1),sNaxes(1),sNaxes(2),OUTDATA,ostat)

666   continue

! close files
      if( ostat == 0 ) then     ! delete output (if any error)
        Call ftclos(26,ostat)
      else
        Call ftdelt(26,ostat) 
      endif
      Call ftclos(25,stat)

      if( (fstat + stat + ostat) /= 0 ) then
        write(*,*) 'The following error(s) has been occured:'
        do
          Call ftgmsg(buf)
          if( buf == ' ' ) exit
          write(*,*) trim(buf)
        enddo
      endif

      nflat = ostat + 1000*stat + 1000000*fstat

    End Subroutine Flat

